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INTRODUCTION 


I_. 

This  report  investigates  the  estimation  of  the  source 
time  parameters  of  underground  nuclear  explosions  from  the 
waveforms  of  short-period  teleseismic  P-waves.  In  the 
simplest  consideration,  and  when  the  source  yield  is  uncon¬ 
strained,  there  are  only  three  source  parameters,  two  that 
describe  the  source  time  function  and  one  for  the  delay  time 
of  the  pP  phase.  There  are,  of  course,  indications  that  the 
pP  arrival  may  not  have  the  same  amplitude  or  shape  as  the 
direct  P  arrival,  presumably  due  to  anelastic  or  non-linear 
effects  between  the  shot  point  and  the  surface.  Thus,  we 
will  also  consider  the  effect  on  the  waveforms  of  a  de¬ 
creased  pP  amplitude.  Another  parameter  which  has  a  signif¬ 
icant  influence  upon  the  observed  waveforms  is  t*  (t*  = 
travel  time/Q  ) .  Accordingly,  the  variation  in  the  source 
parameter  estimates  due  to  a  varying  t*  will  also  be  ad¬ 
dressed. 

There  are  other  phenomena  which  can  influence  the 
waveforms,  such  as  converted  phases  at  the  source  and  re¬ 
ceiver,  of  which  the  latter  seems  to  be  the  most  important. 
We  hope  to  minimize  this  "noise"  principally  by  selecting 
"clean"  stations  when  modeling  the  data.  Figure  1  shows  a 
few  representative  seismograms  of  underground  nuclear 
events,  and  the  characteristic  waveshape  that  is  recorded  by 
short  period  instruments  is  evident.  Also  notice  that  the 
siesmograms  differ  substantially  from  station  to  station 
after  the  first  two  seconds. 


Figure  1.  Representative  WWSSN  short  period  seismograms  for 

an  event  in  Western  Kazakh,  Soviet  Union.  The  time 
scale  is  one  minute  between  time  marks. 
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Before  discussing  the  details  of  determining  the  source 
time  function,  it  is  useful  to  review  the  source  representa¬ 
tion  used.  Based  on  near-field  observations,  Haskell  (1967) 
introduced  a  simple  analytical  formula  for  the  source  dis¬ 
placement  potential  time  history,  with  the  condition  that 
acceleration  at  the  source  is  continuous.  Figure  2  is  taken 
from  Haskell  (1967),  and  shows  the  observations  and  Has¬ 
kell's  analytical  function.  The  far-field  time  function  is 
basically  the  time  derivative  of  the  near-field  potential, 
dropping  the  terms  which  do  not  propagate  to  teleseismic 
distances.  Haskell's  representation  can  be  expressed  in 
analytical  form  as: 

2  3 

Sf-R  =  1  -  e-kt(l+kt4^  +iisn_b(kt)4) 

^  ’  2  6 

There  are  two  parameters  present,  the  time  constant  k  which 
is  proportional  to  the  reciprocal  of  the  far-field  rise 
time,  and  the  overshoot  parameter  b  which  characterizes  the 
amount  of  overshoot,  as  seen  in  Figure  2.  Von  Seggern  and 
Blanford  (1972)  revised  the  Haskell  source  description  by 
allowing  the  near-field  velocity  to  be  discontinuous,  and 
obtained  the  subsequent  analytical  expression: 

=  1  -  ekt(l  +  kt  -  B(kt)2) 

where  once  again  k  is  a  reciprocal  time  constant  and  B 
characterizes  the  overshoot.  As  the  von  Seggern  and  Blan- 


Fig.  2 a.  Granite. 


Tt*C. 

Fig.  2b.  Salt 


Figure  2.  From  Haskell  (1967) 
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ford  form  appears  to  be  consistent  with  teleseismic  spectra, 
this  form  will  be  used  throughout. 

As  indicated  by  Haskell  from  a  simple  scaling  argument, 
-1/3 

k  should  scale  as  Y  '  where  Y  is  the  explosion  yield,  and 
B  is  expected  to  be  fairly  constant  for  a  particular  rock 
unit.  Fitting  the  von  Seggern-Blanford  form  to  the  data  in 
Figure  2,  the  value  of  k  for  the  5  kt.  explosions  varies 
from  ''-8  to  ^16.  Accepting  the  scaling  relation  for  k,  an 
event  of  100  kt.  size  in  those  materials  would  give  k  in  the 
range  ^3  to  ^6.  With  a  more  sophisticated  scaling  law,  von 
Seggern  and  Blanford  used  k  =  16.8  at  Y  =  5  kt.  and  pre¬ 
dicted  that  k  =  9.6  at  Y  =  80  kt. ,  while  the  simple  Y  ~1^3 
scaling  would  give  k  =  6.7  for  Y  =  80  kt. 

As  seen  in  Figure  2,  the  amount  of  overshoot  is  differ¬ 
ent  for  each  event.  The  value  of  B  varies  as:  ~0  for  the 
event  in  tuff,  ~2  for  the  event  in  granite,  and  ~3.3  for  the 
event  in  alluvium.-  It  is  commonly  assumed  that  the  value 
for  B  is  mostly  dependent  upon  the  rock  type  containing  the 
explosion  and  does  not  depend  strongly  on  explosion  yield. 
Placing  an  upper  limit  on  the  B  value  from  knowledge  of  the 
containment  rock  would  be  a  useful  constraint  in  deducing 
the  source  parameters,  as  will  be  seen  later. 

The  delay  time  of  the  surface  reflection  is  dependent 
on  rock  type  and  depth  of  burial.  For  the  common  depths  and 
rock  velocities,  this  time  delay  is  on  the  order  of  a  few 
tenths  of  a  second.  Results  of  non-linear  finite  difference 
calculations,  while  producing  delay  times  greater  than  those 
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predicted  by  elastic  theory,  still  produce  delay  times  on 
the  order  of  tenths  of  seconds  (Mellman  et  al . ,  1980).  As 
teleseismic  short  period  recordings  have  a  duration  of  >1 
sec,  the  pP  arrival  is  entangled  with  the  direct  P  arrival 
and  usually  is  not  visible  as  a  distinct  phase.  Thus,  in 
modeling  the  resultant  waveform  we  need  to  include  the  pP 
contribution.  We  will  assume  the  pP  arrival  has  the  same 
shape  as  direct  P.  This  assumption  is  simply  that  the 
upgoing  and  downgoing  time  functions  are  identical  except 
for  an  amplitude  scale  factor.  Fracturing  of  the  rock  above 
the  explosion  might  cause  the  pP  arrival  to  have  a  broader 
time  function,  thus  delaying  the  peak  amplitude  of  pP 
slightly  beyond  the  time  expected  from  the  burial  depth  and 
rock  velocity.  However,  it  will  be  shown  that  resolving 
different  source  time  paramters  for  pP  is  a  very  difficult 
task. 

It  is  necessary  to  consider  a  range  of  t*  values.  In 
the  earliest  attempts  to  model  teleseismic  seismograms,  it 
was  assumed  that  t*  =  1  sec  and  that  t*  was  constant  for  the 
teleseismic  distance  range  (A^30  to  80  deg.).  This  value 
has  been  in  popular  use  since  then,  if  only  because  there 
was  no  compelling  evidence  for  a  different  value.  More 
recent  work  using  digital  short  period  instruments  peaked  at 
a  higher  frequency  than  the  WWSSN  short  period  (which  is 
peaked  at  1.4  htz )  seems  to  indicate  that  t*  could  be  sub¬ 
stantially  lower  than  1  (e.g.,  Mellman  and  Hart,  1980;  Der, 
et  al.  1979),  perhaps  0.5-0. 6,  and  in  some  cases  values  of 


» 


t 
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0. 1-0.2  have  been  discussed.  It  is  rather  disturbing  to 
have  an  order  of  magnitude  uncertainty  in  a  quantity  that 
appears  in  an  exponential,  as  the  observed  amplitudes  depend 
upon  t*  as, 

A^exp(-rtft* ) 

Thus,  at  f  =  1  hz . ,  allowing  the  range  t*  =  0.1  to  t* 
=  1.0  introduces  a  factor  of  ~17  uncertainty  in  the  ampli¬ 
tude.  As  the  amplitudes  of  the  P-waves  are  used  to  deter¬ 
mine  the  yield,  clearly  the  value  of  t*  is  of  primary  impor¬ 
tance.  These  uncertainties  in  t*  also  affect  the  deter¬ 
mination  of  the  other  source  parameters  as  well. 
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1 1 .  METHOD 

As  a  large  number  of  presumed  underground  nuclear 
events  in  the  Soviet  Union  have  been  well  recorded  on  the 
WWSSN  short  period  instruments,  a  method  has  been  developed 
which  utilizes  these  photographically-recorded  seismograms. 
The  essential  idea  is  that  the  short  period  waveforms  can  be 
entirely  characterized  by  the  peaks  and  troughs  (unless 
there  are  obvious  inflections).  That  is,  the  waveform  can 
be  reproduced  to  within  the  thickness  of  the  photographic 
recording  by  using  a  cosine  interpolation  between  the  peaks 
and  troughs  (see  Figure  1).  The  available  independent  data 
are  the  relative  peak  amplitudes  and  peak  times.  Thus,  for 
the  typical  short  period  waveform  with  four  peaks,  we  can 
represent  the  waveform  with  only  eight  independent  numbers 
corresponding  to  the  amplitudes  and  arrival  times  of  those 
peaks  (see  Figure  3).  This  type  of  representation  has  also 
been  used  by  Somerville,  Wiggins,  and  Ellis  (1976). 

In  previous  studies,  the  time  domain  technique  for 
determining  the  teleseismic  time  function  of  nuclear  events 
was  to  "match"  synthetic  seismograms  with  the  observed 
seismograms  in  a  trial-and-error  fashion  (e.g.,  Burdick  and 
Helmberger,  1979).  Instead  of  applying  this  rather  sub¬ 
jective  method,  we  have  used  a  formal  waveform  inversion 
based  on  the  peak-trough  data  characterization  (note  that 
the  "match"  criteria  is  similar  to  that  used  when  visually 
matching  the  synthetic  to  the  data).  Setting  up  the  formal 
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MODEL 
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k,  B 


lpP 


Source 


Figure  3.  Synthetic  seismograms  are  constructed  by  convolving 
the  source  time  function  with  the  instrument  and 
Q  operator.  A  typical  source  time  function  is 
shown,  along  with  the  instrument  convolved  with  the 
Q  operator  for  t*  =  1.  The  data  and  synthetic  are 
then  parameterized  by  the  peak  amplitudes  and  times 
(as  shown  in  the  lower  portion) . 


t 


10 


inversion  problem  allows  for  an  investigation  of  the  unique¬ 
ness  of  the  solution.  The  trade-off  between  the  three 
parameters  ( k ,  B,  t  )  which  produce  similar  seismograms  has 
not  previously  been  well  explored.  Additionally,  the  method 
could  potentially  be  used  for  the  systematic  estimation  of 
source  parameters. 

In  the  time  domain,  the  seismogram  S(t)  is  constructed 
as  the  convolution, 

S(t)  =  I ( t ) *Q( t; t* ) *F( t ; k, B) *D( t; t  , R) 
where  I  is  the  instrument,  Q  is  the  attenuation  operator,  F 
is  the  source  time  function,  and  D  is  the  half-space  Green's 
function  (a  delta  function  at  t  =  0  followed  by  a  negative 
i  delta  function  at  the  pP  delay  time  t  ),  and  R  is  the  ef- 

fective  surface  reflection  coefficient  which  we  initially 
j  take  to  be  -1.  In  the  teleseismic  range,  the  earth  response 

,  results  only  in  a  spreading  factor,  there  are  no  significant 

i 

waveshape  changes  at  the  periods  of  interest.  In  our  pro¬ 
cedure,  we  fix  I  and  Q,  then  determine  the  function  F*D  as 
parameterized  by  k,  B,  and  t  .  Formally,  as  soon  as  we 
discretize  S(t)  in  the  above  relation  there  is  no  unique 
solution  for  the  unconstrained  function  F*D.  It  is  neces¬ 
sary  to  impose  a  smoothness  condition  on  F*D.  Notice  how- 
!  ever,  that  F  is  parameterized  by  just  two  variables,  k  and 

B,  and,  indeed,  is  itself  a  smooth  function.  Thus,  although 
the  source  parameters  occur  in  a  nonlinear  fashion  in  F,  we 
|  might  hope  that  the  three  source  parameters  could  be  deter- 

i 

[  mined  uniquely. 
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As  the  parameters  appear  in  a  nonlinear  fashion,  we 
adopt  the  usual  linearization  procedure  in  formulating  the 
problem.  Given  initial  estimates  for  the  parameters,  we 


form  the  error  vector  £,  e  =  d  -  S(k,B,tp),  where  d  is  the 

data.  If  the  error  is  unacceptable,  we  want  to  change  S  to 
2 

force  £  toward  zero, 

0  =  d  -  (S  +  6  S) 


then 


6S  =  e 

To  find  t2ie  required  parameter  change,  we  expand  the  syn¬ 
thetic  around  the  current  parameter  values, 


6S 


dS 

dp7Apj+‘ ‘ ' 


Ignoring  the  higher  order  terms,  we  then  have 


dS 


j 


a  linear  system  of  the  form  Ay  =  b,  and  with  three  para¬ 
meters  and  eight  data  values,  an  overdetermined  system. 
After  solving  the  above  system,  the  new  parameter  values 
will  be  ^  p  =  q£  +  Ap.  As  the  higher  terms  in  <5S  have  been 
ignored,  and  data  errors  lead  to  an  incompatible  system,  it 
is  likely  that  the  error  has  not  been  reduced  completely  to 
zero,  and  =  d  -  S(^p).  Thus,  we  iterate  with  the  same 
procedure  to  reduce  n£  to  the  desired  level,  dictated  by  the 
data  variance.  In  most  inversion  procedures  where  there  can 
be  large  numbers  of  parameters  and  data,  one  would  select 


» 
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different  starting  models  (Qp),  and  allow  the  method  to 
iterate  to  test  whether  the  procedure  converges  to  the  same 
solution.  In  our  situation,  given  the  small  size  of  the 
system  (8x3),  it  is  feasible  to  "map  out"  the  three  dimen¬ 
sional  parameter  space,  thereby  examining  the  properties  of 
convergence  for  synthetic  data.  One  display  method  is  to 
contour  the  modulus  squared  of  the  error  vector  as  a  func¬ 
tion  of  the  parameters, 

2  2 
£(£l  =  d-s(p) 

However,  contouring  the  length  squared  error  vector  does  not 
display  the  length  of  the  parameter  perturbation,  particu¬ 
larly  when  damping  is  used.  Hence,  the  properties  of  con¬ 
vergence  are  displayed  here  with  "vectorgrams" .  A  vector- 
gram  plot  the  parameter  perturbations  as  a  function  of  the 
parameter  values.  At  a  particular  p,  let  A^  j  =  dS^/dp^ , 
then  solving  the  system  A  —  Apj  =  we  obtain  Ap^  =  A_1ji 
£^.  Thus  at  any  p  we  plot  the  associated  Ap,  with  the  tail 
of  the  vector  at  p  and  the  head  at  p  +  Ap.  We  can  then  plot 
the  Ap  at  network  of  values  for  p.  Figure  4  shows  a  vec- 
torgram  using  synthetic  data  (in  which  the  "true"  solution 
is  of  course  known).  The  three-dimensional  network  of 
parameter  values  samples  k  and  B  quite  well,  covering  the 
range  of  1  to  9  and  0  to  4  respectively,  while  only  three 
values  of  t  (0.4,  0.5,  0.6  sec.)  are  sampled.  The  limited 
sampling  of  t  is  adequate  as  At  is  well  resolved  every- 

r  r 

where,  which  is  demonstrated  in  a  later  section. 
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Figure  4 
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\fectorgram  for  synthetic  data.  The  data  are  taken  from  a 
synthetic  seismogram  computed  with  the  following  parameters 
k=5,  B=2,  tp=0.5,  and  t*=l.  The  vectorgram  displays  the 
parameter  perturbations  calculated  at  a  network  of  (k,B,tp) 
values,  with  1%  damping.  Starting  at  initial  values  for 
k,  B,  and  t  ,  the  vectorgram  shows  the  path  that  the  inver¬ 
sion  procedure  would  follow. 
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III.  RESULTS 

3 . 1  Properties  of  the  Inversion  Method 

There  are  several  features  in  Figure  4  that  persist  for 

I 

many  of  the  cases  that  have  been  considered.  One  important 
feature  in  Figure  4  is  that  the  method  does  converge  to  the 
"true"  solution,  if  a  suitable  starting  place  is  selected. 
If  one  starts  with  small  values  of  k  and  B,  the  method 
proceeds  directly  toward  the  answer.  Figure  5  is  a  sample 
inversion  run  starting  with  small  values  of  k  and  B.  How¬ 
ever,  when  starting  at  any  other  corner  of  the  k,B  plane, 
the  method  would  proceed  to  the  "trough"  and  stop  without 
reaching  the  "true"  answer.  This  trade-off  curve,  which 
extends  toward  a  larger  k  and  a  larger  B  from  the  correct 
solution,  appears  to  represent  an  inherent  nonuniqueness  in 
the  solution,  and  is  present  in  all  vectorgrams  with  artifi¬ 
cial  data  for  a  variety  of  short-period  instruments.  To 
show  that  this  trough  is  not  an  artifact  of  the  inversion 
procedure,  Figure  6  compares  the  synthetic  seismogram  for 
parameters  (5.0,  2.0,  0.5)  with  two  other  synthetic  seismo¬ 
grams  computed  with  parameter  values  at  either  end  of  the 
trade-off  curve,  and  indeed  they  are  quite  similar.  Looking 
at  the  source  description  formula,  the  nature  of  this  trade¬ 
off  is  not  obvious.  If  we  compare  the  spectral  amplitudes 
of  the  source  time  functions  (Figure  7),  it  would  seem  that 
the  k,B  trade-off  combinations  maintain  the  level  of  the 
high  frequency  slope  relative  to  the  low  frequency  level. 

It  appears  that  the  height  of  the  spectral  high  due  to  the 


15 


Figure 

» 


STARTING 

MODEL 


k  =  2.00 


B  -  0.50 

tpp  ~  0. 60 


FIRST 

ITERATION 


k  -  2.94 
B  -  1.35 


*pP 


=  0.53 


SECOND 

ITERATION 


k  -  3.99 
B  -  1.83 


*pP 


-  0.51 


V 


THIRD 

ITERATION 


k  *  4.48 
8  *  2.08 
tpp  -  0.51 


F 


5  sec 


5.  An  inversion  sequence  for  the  synthetic  data  used  in  con¬ 
structing  the  vectorgram  in  Figure  4 .  The  source  time  func-  I 

tions  and  resultant  synthetic  seismograms  are  shown  for  i 

the  starting  model  and  three  iterations.  The  solid  trace  at 
each  iteration  is  the  "data"  (computed  with  k=5,  B-2 ,  t  =0.5), 
and  the  dashed  trace  is  the  synthetic  corresponding  to  the  [ 

source  parameters  of  that  iteration.  There  is  an  excellent  j 

agreement  after  three  iterations.  j 

i 

.  ..  | 


Comparison  of  synthetics  along  the  trade-off  curve, 
using  the  WWSSNSP  instrument  with  t*=1.0  and  tpp=0. 
sec  in  all  cases.  The  synthetic  seismogram  for  k=5 
and  B=2  is  shown  in  both  (a)  and  (b)  as  the  solid 
trace.  The  dashed  trace  in  (a)  is  the  synthetic 
seismogram  for  k=4.4  and  B=3.5.  The  dashed  trace 
in  (b)  is  the  synthetic  seismogram  for  k=6.5  and 
B=1 . 3 . 


SPECTRAL  AMPLITUDE 


k  -  5.0,  B  =2.0 


k  =  4.4,  B  =  3.5 


Figure  7.  The  spectral  amplitudes  for  three  different  (k,B)  combina¬ 
tions.  These  three  combinations  are  on  the  trade-off 
curve  in  Figure  4.  It  appears  that  k,B  values  on  the  trade 
off  curve  are  characterized  by  maintaining  the  high  fre¬ 
quency  decay  slope  at  the  same  level  relative  to  the  zero 
frequency  level. 


overshoot  does  not  exert  a  strong  influence  on  the  wave¬ 
forms.  This  may  be  partly  due  to  the  spectral  null  that  is 
introduced  at  about  this  frequency  by  the  pP  arrival. 


Based  on  other  examples  of  synthetic  data,  we  conclude 
that  this  trade-off  curve  is  a  persistent  feature  of  this 
inversion  technique,  but  that  the  "true"  solution  lies  at 
the  corner  of  the  trade-off  curve.  Thus,  starting  with  low 
values  of  k  and  B,  then  allowing  the  inversion  technique  to 
proceed  should  provide  at  least  systematic  estimates  of  the 
source  parameters. 

It  was  stated  earlier  that  t  is  well  resolved,  in  that 
the  correct  value  of  t  is  closely  approximated  after  just 
one  step  from  a  large  part  of  the  k,B,tp  network.  Once 
again,  many  different  vectorgrams  for  different  cases  exhib¬ 
it  this  feature.  In  this  regard,  it  is  instructive  to  plot 
the  resolution  matrices  at  each  gridpoint.  The  resolution 
matrix,  R,  is 

R  =  A_1A 


where  A-1  is  the  particular  inverse  used.  In  the  case  of  no 
damping  and  no  data  weighting,  the  LANCZOS  (1961)  inverse 
is, 


.-1 


=  (ATA)"1AT 


With  this  definition  R  is  the  identity  matrix,  that  is,  all 
three  parameters  are  resolved,  the  three  diagonal  elements 
equaling  1.  However,  the  solutions  are  somewhat  unstable. 
Moving  to  an  adjacent  gridpoint  results  in  a  significantly 
different  Ajd.  It  was  found  that  a  1%  damping  gave  the 
desired  stability,  resulting  in  the  inverse, 


l 
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A"1  =  [ATA  +  cTr(ATA)]“1AT 

With  this  form  of  A  ',  R  is  no  longer  the  identity  matrix, 
and  the  parameters  associated  with  a  small  eigenvalue  will 
have  a  diagonal  element  between  0  and  1.  Ti  three  diagonal 
elements;  corresponding  to  three  parameters  k,B,  and  t^;  are 
plotted  as  histograms  at  each  gridpoint  (Figure  8),  with  the 
top  of  each  box  corresponding  to  a  value  of  1.  It  is  quite 
obvious  that  At^  is  resolved  everywhere,  and  this  is  typical 
for  other  t  planes.  There  are  two  other  important  conclu¬ 
sions  to  be  reached  from  Figure  8.  One  is  that  the  value  of 
B  is  only  well  resolved  for  small  values  of  B,  that  is,  the 
waveform  "knows"  whether  B  =  0  or  B  >0,  but  beyond  a  certain 
limit  the  value  of  B  is  not  well  determined.  Also,  the 
value  of  k  is  well  resolved  for  small  k,  up  to  k-^5.  Beyond 
that,  the  value  of  k  is  poorly  determined.  This  can  be 
understood  as  for  k>5,  the  basic  source  function  pulse  width 
is  <0.5  sec.,  approximately  the  minimum  value  of  any  peak- 
to-peak  time.  Thus,  in  the  limit  of  large  k,  the  source 
function  is  basically  a  delta  function  compared  to  the  tele- 
seismic  waveform  (i.e.,  there  would  be  negligible  differ¬ 
ences  in  the  waveforms  between  k  =  20  and  k  =  30).  This  is 
an  inherent  difficulty  in  characterizing  the  data  by  the 
peaks,  very  little  information  can  be  ascertained  for  per¬ 
iods  substantially  smaller  than  the  period  associated  with 
the  minimum  peak  time  value.  It  should  be  pointed  out  that 
the  quantitative  values  of  these  limitations  is  a  function 
of  the  instrument  response.  The  values  discussed 
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Figure  8.  The  diagonal  elements  of  the  resolution  matrices  for  the  network  of  (k,B) 

values,  with  t  =0.5.  At  each  (k,B)  point,  a  histogram  displays  the  size  of 
the  three  diagbnal  elements,  the  first  bar  is  for  k,  the  second  bar  is  for 
and  the  third  bar  is  for  tp.  The  value  at  the  top  of  each  histogram  is  1, 
i.e.,  well  resolved,  and  the  value  at  the  bottom  is  0. 
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above  are  those  appropriate  for  the  WWSSN  short-period 
seismometer.  A  finer  time  resolution  would  be  obtained  from 
digital  recordings  of  higher  frequency  instruments. 

3 . 2  Source  errors  due  to  improper  t*  and  pp  amplitude. 

The  vectorgram  in  Figure  4  is  computed  with  t*  =  1.0, 
as  is  the  synthetic  used  for  data.  Given  the  uncertainty  in 
t*,  we  would  like  to  know  what  are  the  errors  in  the  source 
parameters  associated  with  an  improper  t*  value.  Figure  9 
shows  an  example  where  the  synthetic  data  are  for  a  t*  = 
0.8,  but  we  have  assumed  t*  =  1.0  in  constructing  the  vec¬ 
torgram.  The  form  of  the  vectorgram  is  the  same  as  Figure 
4,  though  the  "true"  solution  (at  corner  of  the  trade-off 
curve)  has  been  shifted  toward  larger  values  of  k  and  B,  and 
tp  is  modified  slightly.  Figure  10  is  an  example  where  the 
synthetic  data  are  for  t*  =  1.2,  and  the  vectorgram  is 
constructed  with  t*  =  1.0.  The  solution  is  then  shifted 
toward  lower  values  of  k  and  B.  Figure  11  shows  the  syn¬ 
thetic  seismograms  for  different  t* ' s  and  source  parameter 
combinations,  demonstrating  that  the  inversion  procedure 
does  find  a  satisfactory  match  for  quite  a  range  of  mismatch 
in  t*,  over  the  k,B  values  considered. 

Figure  12  shows  the  vectorgram  for  synthetic  data 
computed  with  the  proper  t*,  but  with  the  pP  amplitude  =  -.8 
while  the  vectorgram  is  constructed  for  pP  amplitude  =  -1. 

A  shift  similar  to  that  due  to  an  incorrect  t*  is  intro¬ 
duced.  This  may  not  be  a  severe  problem  as  the  theoretical 


22 


Figure  9.  This  vectorgram  shows  the  error  in  parameter  estimates 
introduced  when  the  assumed  t*  is  incorrect.  The  syn¬ 
thetic  data  is  computed  for  k=5,  B=2,  tp=0.5  and  t*=0.8. 
The  vectorgram  is  constructed  assuming  t*=1.0.  Notice 
that  the  "true"  solution  has  been  shifted  from  k=5  and 
B=2 . 
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Figure  10.  This  vectorgram  shows  the  shift  in  parameter  estimates 
with  an  incorrect  t*  value.  The  synthetic  data  is  com¬ 
puted  for  k=5 ,  B=2 ,  tp~0 . 5 ,  and  t*=1.2,  while  the  vector- 
gram  is  constructed  assuming  t*=1.0. 


Figure  11.  Comparison  of  synthetic  seismograms  with  different 

values  for  t* .  In  (a),  the  solid  trace  is  computed  for 
t*=0 . 8 ,  k=5,  B=2,  and  tpp=0.5,  while  the  dashed  trace 
is  for  t*=1.0,  k=6.37,,  B=2.3,  tpP=0.45.  In  (b) ,  the 
solid  trace  is  computed  for  t*=0.6,  k=5,  B=2,  and  tpp=0 
and  the  dashed  trace  is  for  t*=1.0,  k=9.5,  B-4.5.  and 
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Figure  12.  This  vectorgram  shows  the  shift  in  parameter  estimates 
when  the  pP  amplitude  is  incorrect.  The  synthetic 
data  is  computed  with  k=5,  B=2,  tp=0.5,  t*=1.0,  and  pP 
ampli  tude=-0 . 8 ,  while  the  vectorgram  is  computed  for 
pP  amplitude--! . 0  (the  correct  t*  value  is  assumed) . 
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modeling  indicates  that  a  pP  amplitude  £-.85  is  adequate  for 
many  cases  and  this  value  could  be  used  as  the  standard  pro¬ 
cedure.  Regarding  the  possibility  of  allowing  a  different 
source  function  for  pP,  it  seems  clear  from  the  resolution 
discussion  that  these  parameters  would  not  be  resolved  in 
general . 

3 . 3  Results  for  Western  Kazakh  Test  Site 

Thus  far,  we  have  considered  only  error  free  data.  In 
using  real  data  there  are  two  basic  sources  of  noise;  in¬ 
trinsic  record  noise  and  waveform  changes  due  to  converted 
phases  as  discussed  earlier.  We  can  improve  the  signal  to 
record  noise  ratio  by  selecting  larger  events  to  model. 
However,  the  converted  phases  (frequently  referred  to  col¬ 
lectively  as  receiver  structure)  will  persist.  To  gain  some 
insight  into  the  waveform  scatter,  three  events  (Dec.  6, 
1969;  Dec.  12,  1970;  Dec.  23,  1970)  which  occurred  in  the 
Western  Kazakh  region  were  examined.  These  presumed  ex¬ 
plosions  were  located  quite  close  to  each  other  and  the 
estimated  yields  range  from  100  to  240  kt.  (Dahlman  and 
Israelson,  1977).  When  using  all  of  the  North  American 
WWSSN  stations  the  total  scatter  in  the  waveforms  is  quite 
large,  though  the  peak  times  are  in  fair  agreement,  the 
scatter  in  peak  amplitudes  is  distressing,  particularly  the 
fourth  peak.  Directly  comparing  the  waveforms  quickly 
isolates  the  clean  or  relatively  transparent  stations.  The 
data  vector  is  then  formed  by  averaging  together  the  peak 
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values  from  only  a  few  of  the  best  stations.  Figure  13 
shows  two  vectorgrams  for  two  different  averages,  one  for 
the  event  of  Dec.  12,  1970  and  another  for  the  combined 

average.  A  t*  of  1  was  used  in  both  cases.  In  the  top 
vectorgram,  the  data  produces  a  pattern  quite  similar  to 
those  of  the  synthetic  data  discussed  earlier.  The  tradeoff 
curve  is  apparent  and  we  could  immediately  choose  the  best 
solution  at  the  corner  of  the  trade-off  curve  (k~5.5,  B-v2.0, 
t  =  0.5).  The  data  vector  used  for  Kazakh  average  is  quite 
similar  to  that  for  Kazakh  12-12-70,  the  most  significant 
difference  being  the  peak  value  for  the  fourth  peak,  and  yet 
the  trade-off  structure  has  disappeared  and  the  vectors 
point  off  to  the  upper  right-hand  corner.  In  this  situation 
it  is  important  to  consider  the  data  errors.  The  contours 
in  Figure  13  refer  to  the  modulus  of  the  error  vector.  In 
these  units,  a  scatter  in  the  peak  times  of  0.05  sec  and  of 
^10%  in  the  relative  peak  amplitudes  gives  an  error  modulus 
of  ^0.4,  and  we  cannot  expect  the  inversion  to  reduce  the 
error  much  below  this  level.  The  closed  0.2  contour  for 
Kazakh  12-12-70  requires  the  peak  times  to  be  within  ~.03 
sec.  and  the  amplitudes  to  be  matched  to  ~5%.  Although 
there  is  no  closed  contour  for  Kazakh  average,  it  would  be 
reasonable  to  stop  the  inversion  procedure  at  the  0.4  con¬ 
tour.  In  this  case,  it  is  of  no  great  concern  that  there  is 
no  closed  contour  as  this  is  probably  due  to  the  data  scat¬ 
ter  causing  incompatibility.  In  this  regard,  the  presence 
of  the  trade-off  curve  for  Kazakh  12-12-70  might  be  fortui- 
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tous  as  opposed  to  better  data  quality.  Figure  14  shows  an 
inversion  run  using  the  Kazakh  12-12-70  data,  and  clearly  a 
good  fit  is  obtained  after  three  iterations,  though  it  has 
not  yet  progressed  to  the  trade-off  curve  corner. 

To  quickly  consider  other  t*  values,  Figure  15  plots  a 
vectorgram  using  the  Kazakh  12-12-70  data  and  a  t*  =  0.6. 
As  expected,  the  source  parameters  move  to  a  smaller  k  and 
slightly  smaller  B,  the  trade-off  curve  is  present,  and  the 
mismatch  acceptable.  Lowering  the  t*  moves  the  source 
parameters  into  the  region  of  resolution.  Also  note  that  t 
appears  to  be  well  resolved  in  these  cases. 

3.4  Inversion  using  SRO  short  period  instruments 

To  extend  the  resolution  of  the  source  time  parameters 
it  is  necessary  to  use  seismograms  with  a  higher  frequency 
content.  The  SRO  short  period  instrument  is  peaked  at  ~2.5 
hz.,  in  contrast  to  the  ~1.3  hz.  peak  of  the  WWSSN  instru¬ 
ment.  Therefore,  we  might  expect  to  increase  the  resolution 
with  the  SRO  instrument.  Also,  the  SRO  is  digitally  re¬ 
corded  producing  less  ambiguity  in  determining  the  peaks 
(indeed  ultimately  we  may  wish  to  use  more  than  just  the 
peaks).  To  demonstrate  that  the  basic  properties  of  the 
inversion  method  apply  when  using  the  SRO  instrument,  Figure 
16  shows  the  vectorgram  for  synthetic  data  with  t*  =  1.0. 
Notice  that  t  is  slightly  less  stable  than  for  the  WWSSN 
instrument.  Figure  17  plots  the  diagonal  elements  of  the 
resolution  matrices  for  t*  =  1.0.  It  is  somewhat  discour- 
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Figure  14.  An  inversion  sequence  for  the  Western  Kazakh  event  of 
12-12-70.  As  in  Figure  5,  the  source  time  functions 
are  shown  on  the  left  and  the  data  (solid  trace)  and 
synthetic  (dashed  trace)  are  shown  on  the  right.  An 
acceptable  match  is  obtained  after  three  iterations. 
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Figure  17.  The  diagonal  elements  of  the  resolution  matrices  for  the  SRO  short  period 

instrument,  computed  for  t*=1.0.  As  in  figure  8,  each  histogram  depicts  the 
resolution  of  the  three  parameters;  k,  B,  and  t_,  ranging  from  0  to  1 . 
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aging  that  we  have  not  increased  the  resolution  substantial¬ 
ly,  but  this  is  verified  as  the  SRO  synthetics  for  t*  =  1.0 
have  peak  times  similar  to  the  WWSSN  synthetics.  Clearly, 
the  broad  attenuation  operator  for  t*  =  1.0  is  limiting  the 
resolution. 

SRO  data  at  several  stations  has  been  examined,  and 
CTAO  was  chosen  as  one  of  the  clean  stations  recording 
events  from  the  Semipalatinsk  area.  When  attempting  to 
model  the  CTAO  seismograms,  it  was  found  that  t*  =  1.0  was 
not  acceptable,  as  the  peak  to  peak  times  of  the  synthetic 
seismograms  were  too  long,  even  as  k  became  quite  large  and 
t  became  quite  small.  Thus,  one  is  compelled  to  use  a 
lower  t*  to  produce  synthetics  with  the  appropriate  pulse 
widths.  Though  it  is  difficult  to  determine  exactly  what 
value  of  t*  is  appropriate,  t*  =  0.4  appeared  to  be  adequate 
for  the  seismograms  considered.  Another  problem  encountered 
was  the  value  of  the  pP  amplitude.  Testing  various  values, 
the  most  compatible  amplitude  is  -1.3.  This  is  a  bit  dif¬ 
ficult  to  reconcile,  and  the  present  interpretation  is  that 
this  value  includes  a  converted  phase. 

We  have  divided  the  events  recorded  by  CTAO  into  two 
categories  based  upon  size.  Recalling  the  earlier  dis¬ 
cussion  on  source  parameter  scaling,  it  would  be  interesting 
if  we  could  resolve  such  differences  at  an  individual  sta¬ 
tion.  The  events  that  will  be  considered  are:  6-23-79 
relative  amplitude  =  246,  6-11-78  relative  amplitude  =  128, 
and  8-4-79  relative  amplitude  =  215.  As  6-23-79  is  only  a 
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factor  of  2  larger  than  6-11-78,  we  would  expect  that  the 
value  of  k  for  6-23-79  will  be  smaller  by  only  ^20%.  As  the 
distortion  in  the  fourth  peak  is  too  large  to  ignore,  the 
fourth  peak  amplitude  and  time  have  been  given  the  weight  of 
0.01.  Figure  18  shows  the  k,B  paths  followed  by  the  dif¬ 
ferent  events.  It  appears  that  the  smaller  6-11-78  event 
consistently  prefers  a  slightly  larger  value  for  k  and  B  as 
the  iterations  proceed.  For  example,  starting  at  (2.0,  .25, 
.5);  after  two  iterations  the  6-23-79  k  value  is  ^12% 
smaller  than  the  6-11-78  value,  after  four  iterations  is  is 
~14%  smaller,  and  after  five  iterations  it  is  ^15%  smaller. 
These  values  are  roughly  consistent  with  the  expected  dif¬ 
ference  based  on  the  simple  scaling  argument,  though  we 
cannot  place  a  strong  emphasis  on  this  result  until  we 
understand  the  causes  of  the  mismatch  at  CTAO. 

Perhaps  of  more  interest  is  tp .  Figure  19  plots  the  tp 
values  for  the  three  events,  and  despite  the  oscillatory 
behavior  there  appears  to  be  a  resolvable  difference  in  the 
pP  time  for  6-11-78,  in  the  direction  consistent  with  a 
shallower  depth  of  burial  for  the  smaller  event.  Figure  20 
shows  the  inversion  results,  with  the  final  time  functions 
and  synthetics  for  the  two  events;  6-23-79  and  6-11-78, 
recorded  at  CTAO.  The  differences  in  the  seismograms  are 
subtle,  but  can  be  seen  with  the  aid  of  a  ruler. 

Concluding  this  preliminary  investigation  of  SRO  seis¬ 
mograms,  some  of  the  source  to  receiver  combinations  require 
a  t*  <1  to  match  the  peak  to  peak  times,  and  using  the 


The  (k,B)  paths  followed  by  the  inversion  procedure.  The  data  used  are 
three  events,  6-23-79,  6-11-78,  and  8-4-78,  recorded  by  the  SRO  short  period 
station,  STAO.  The  t*  value  is  0.4.  Different  starting  models  are  used. 

The  parameter  estimates  for  6-11-78  are  consistently  different  than  for 
6-23-79  and  8-4-78. 
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(a) 

INVERSION  RESULT  FOR  6-23-79,  CTAO 


(b) 

INVERSION  RESULT  FOR  6-11-78,  CTAO 


sec 


Figure  20.  Inversion  results  for  the  two  events:  (a)  6-23-79  and 
I  (b)  6-11-78,  recorded  by  the  CTAO  SRO  short  period  instru¬ 

ment  (the  solid  trace  in  (a)  and  (b) .  With  initial  values 
of  k=2,  B=2 . 5 ,  and  tpp=0.5  sec,  the  time  functions  and 
synthetics  (dashed  trace)  after  five  iterations  are  shown. 
The  source  values  for  (a)  are  k=5.1,  B=3.1,  and  tpp=.41. 

A  t*=0.4  and  surface  reflect  coeff=-1.3  were  used. 
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records  at  CTAO,  there  is  a  suggestion  that  the  relative 
value  of  k  between  different  events  might  be  resolvable  to 
<20%,  and  the  value  of  t  might  be  resolvable  to  <.05  sec. 

I  ** 
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